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Abstract 

We present a parameter study of simulations of fragmentation regulated by gravity, magnetic 
fields, ambipolar diffusion, and nonlinear flows. The thin-sheet approximation is employed with 
periodic lateral boundary conditions, and the nonlinear flow field ("turbulence") is allowed to 
freely decay. In agreement with previous results in the literature, our results show that the onset 
of runaway collapse (formation of the first star) in subcritical clouds is significantly accelerated 
by nonlinear flows in which a large-scale wave mode dominates the power spectrum. In addition, 
we find that a power spectrum with equal energy on all scales also accelerates collapse, but by a 
lesser amount. For a highly super- Alfvenic initial velocity field with most power on the largest 
scales, the collapse occurs promptly during the initial compression wave. However, for trans- 
Alfvenic perturbations, a subcritical magnetic field causes a rebound from the initial compres- 
sion, and the system undergoes several oscillations before runaway collapse occurs. Models that 
undergo prompt collapse have highly supersonic infall motions at the core boundaries. Cores in 
magnetically subcritical models with trans- Alfvenic initial perturbations also pick up significant 
systematic speeds by inheriting motions associated with magnetically-driven oscillations. Core 
mass distributions are much broader than in models with small-amplitude initial perturbations, 
although the disturbed structure of cores that form due to nonlinear fiows does not guarantee 
subsequent monolithic collapse. Our simulations also demonstrate that significant power can (if 
present initially) be maintained with negligible dissipation in large-scale compressive modes of 
a magnetic thin sheet, in the limit of perfect flux freezing. 
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1. Introduction 



Magnetic fields and supersonic turbulence are two mechanisms that are commonly in- 
voked for the regulation of star formation in o ur Galaxy to the observationally estimated 
rate of - 3 - SM© yr"! (see iMcKeel . Il989t ). This is at least one hundred times less 
than the rate implied by the gravitational fragmentation timescale of the molecular gas 
in the Galaxy calculated from its mean density. Put another way, the global Galactic 
star formation efficiency (SFE) is about 1% of molecular gas per free-fall time. Interest- 
ingly, the same SFE als o applies to nearby ind ividual star-forming regions such as the 



Taurus molecular cloud ( Goldsmith et al. The relatively low Galactic SFE is one 



fundamental constraint on the global properties of star formation. The existence of a 
broad-tailed core mass function (CMF) that is a lognormal with a possible power-law 



high-mass tail, is another. In fact, the observed form of the CMF (e.g. lMotte et al.L 1X9981 ) 



is similar to that of the stellar initial mass function, the IMF. Other important star forma- 
tion cons traints specific ally appl ying to cores include the generally subsonic relative infa ll 
motions fTa falla et all ., 1998: Williams et "all. 1 19991: [T ee et al.l.l200lHCasem et"alll2002l) 
the l ow (subsonic or transonic) systematic core speeds ( Andre et all 2007t Kirk et al 
2007 ). and the somewhat non-circular projected shapes ( Mvers et al.l . 199X1 : Jones et al 



20011 ) . The relatively low speeds are an important property since the cores are embedded 



in molecular clouds whose overall internal random motions are highly supersonic. 

Since most stars form in clusters or loose groups, it seems clear that some sort of 
fragmentation process is at work in those regions. There are several qualitatively distinct 
modes of fragmentation to consider. The simplest process is gravitational fragmentation, 
which can be divided into cases with mass-to-flux ratios that are supercritical (gravity- 
dominated) and subcritical (ambipolar-diffusion-driven). An alternate and distinct star 
formation mechanism is turbulent fragmentation, dominated by nonlinear flows, which 
can also occur in clouds with supercritical and subcritical mass-to-flux ratios. The limit 
of highly supercritical clouds also corresponds to super- Alfvenic turbulence in the case 
that turbulent an d gravitational energi e s hav e comparable magnitude. This limit has 
been advocated bv IPadoan fc Nordlundl (|2002l) . However, we favor the transcritical and 
trans-Alfvenic cases on theoretical and observational grounds, as discussed in Section 
4. For a more complete understanding of star formation, a study of the interplay of 
magnetic fields, turbulence, and ambipolar diffusion is therefore of great importance. In 
this paper, we carry out an extensive parameter survey of magnetic field strengths and 
initial nonlinear perturbations, facilitated by our use of the thin-sheet approximation, 
and discuss the consequences of our results. Such broad parame ter studies remain out 
of reach for fully thre e-dimensional non-ideal MHD simulations ( Kudoh fc Basu . 2008t 
Nakamura fc"Lil . l2008h . 



In a previous paper (|Basu et al.ll2009L hereafter BCW; see also iBasu fc Ciolekll2004l) . 
we studied the nonlinear evolution of gravitational fragmentation initiated by small- 
amplitude perturbations, including the effects of magnetic fields and ambipolar diffusion. 
An extensive parameter study was performed, encompassing the supercritical, transcrit- 
ical, and subcritical cases, and also accounting for varying levels of cloud ionization and 
external pressure. Some main findings of that paper were th at fragment spacings in the 
nonlinear phase agree with the predictions of linear theory ( Ciolek fc Basu . 2006I ). and 
that the time to runaway collapse from small-amplitude white-noise initial perturbations 
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is up to ten times the growth time calculated from linear theory. ICiolek fc Basu (|2006l) 
showed that transcritical gravitational fragmentation can lead to significantly larger size 
and mass scales than either the supercritical or subcritical limits. BCW found that CMFs 
for regions with a single uniform initial mass-to- flux ratio are sharply peaked, but that the 
sum of results from simulations with a variety of initial mass-to-flux ratios near the crit- 
ical value can produce a broad distribution. This represents a way to get a broad CMF, 
of the type observed, without the need for nonlinear forcing. Importantly, only a narrow 
initial distribution of initial mass-to-flux ratio is needed to create a relatively broad CMF. 
Additionally, BCW showed that different ambient mass-to-flux ratios in different regions 
lead to observationally distinguishable value s of infall motions, core shapes, and mag- 
netic field line curvature. iKudoh et al. ( 2007 ) have recently performed three-dimensional 
simulations of gravitational fragmentation with magnetic fields and ambipolar diffusion, 
and verified some of the main findings of the thin-sheet models. Particularly, they also 
found the dichotomy between subsonic maximum infall speeds in subcritical clouds and 
somewhat supersonic speeds in supercritical clouds. 

The inclusion of nonlinear (hereafter, "turbulent") initi al conditions to fragrn enta- 
tion models incl u ding ambipolar diffusion was introduced bv lLi fc Nakamural (|2004r ) and 
Nakamura fc Li ( 20051 ). employing the thin-sheet approximation. They found that the 



timescale of star for mation was reduced sig nificantly by the initial motions with power 
spectrum v'^ (x fc"'* ( Li fc Nakamural 2004), to become ^ 10^ yr for an initially some- 
what subcritical cloud rather than ~ 10^ yr. By continuing to integrate past the collapse 
of the first core through the use of an artificially stiff equation of state for surface densi- 
ties 10 times greater than the initial value, they found that magnetic fields nevertheless 
prevented most material from collapsing to form stars. iKudoh fc Basu (2008) have veri- 
fied that the mode of turbulence-accelerated magnetically-regulated star formation also 
occurs in a fully three-dimensional simulation. While three-dimensional simulations are 
resour ce-limited, and large parameter studies cannot yet be performed, Kudoh fc Basul 
(j2008l ) showed that this mode of star formation proceeds though an initial phase of en- 
hanced ambipolar diffusion created by the small length scale generated by the large-scale 
compression associated with the initial perturbation. If this is not sufficient to raise the 
maximum mass-to-flux ratio above the critical value, there is a rebound to lower den- 
sities. However, the highest density regions remain well above the initial mean density, 
and here the ambipolar diffusion proceeds in a quasis tatic manner (ex Pn assum ing 



force balance between gravity and magnetic forces - see Mouschovias fc Ciolek 1999f) but 
at an enhanced rate due to the raised density. 

In this paper, we study the effect of large-amplitude nonlinear initial perturbations on 
the evolution of a thin sheet whose evolution is regulated by magnetic fields and ambipolar 
diffusion. We focus on the early stages of prestellar core formation and evolution, and 
do not integrate past the runaway collapse of the flrst core. Therefore, the effects of 
protostellar feedback through outflows do not need to be added. These simplifications 
allow us to run a large number of simulations. We perform an extensive parameter study 
and also study many realizations of models with a single set of parameters, since the initial 
turbulent state is inherently random. Some important questions that we can address and 
which have not been answered in previous papers are as follows. In which parameter space 
do nonlinear velocity fields lead to prompt collapse, and in which cases can magnetic fields 
cause a rebound from the first compressions? How sensitively does the time until runaway 
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collapse depend upon the values of different parameters? What is the effect of different 
power spectra of perturbations? Is there a qualitative difference between Alfvenic and 
super- Alfvenic perturbations? How do velocity profiles in the vicinity of cores vary in the 
different scenarios? What are the systematic speeds of cores? 

Our paper is organized in the following manner. The model is described in Section 
2, results are given in Section 3, and a discussion of results, including speculation and 
implications for global star formation in a molecular cloud are given in Section 4. We 
summarize our results in Section 5. 



2. Physical and Numerical Model 

We employ th e magnetic thin-sheet approximation, as laid out i n detail in severa l 
previous papers (ICiolek fc MouschoviasL IiQqI ICiolek fc Basul . 120061 : iBasu et all . [200i). 



Physically, we are modeling the dense mid-layer of a molecular cloud, and ignoring the 
more rarefied envelope of the cloud. We assume isothermality at all times. The basic 
equations governing the evolution of a model cloud (conservation of mass and momentum. 
Maxwell's equations, etc.) are integrated along the vertical axis from z = —Z{x,y) to 
z = +Z(x,y). In doing so, a "one-zone approximation" is used, in which all quantities 
are taken to be independent of height within the sheet. The volume density of neutrals 
Pn is calculated from the vertical pressure balance equation 

Pncl = ^Gal^ Pext + \ " , (1) 

where Cg is the isothermal sound speed, a-a{x^ y) = /_^^ Pn{x, y) dz is the column density 
of neutrals, Pext is the external pressure on the sheet, and and By represent values 
of magnetic field components at the top surface of the sheet, z = +Z. 

We solve normalized versions of the magnetic thin-sheet equations. The unit of velocity 
is taken to be Cs, the column density unit is (Jn.Q, and the unit of acceleration is 27rGCTn,o, 
equal to the magnitude of vertical acceleration above the sheet. Therefore, the time 
unit is to = Cs/27rGan,o, and the length unit is Lq — Cg/27rGcrn,o- From this system 
we can also construct a unit of magnetic field strength, Bq = 27rG"'^/^(7n.o- The unit of 
mass is Mq — cf/(47r^G^ crn,o)- Here, (Jn,o is the uniform neutral column density of the 
background state. Typical values of our units are given in the Appendix. With these 
normalizations, the equations used to determine the evolution of a model cloud are 

-Vp • ((Tn ■»„) , (2) 
-Vp • (o-nt^nVn) + Ft + Fm + CFnOp, (3) 

-yp-{B,,eqVi), (4) 

-CoffVpCTn, (5) 

B,,eq {Bp - Z VpB,,,^) + 0{VpZ), (6) 

T'ni^O / Pn.oV' IP l'7\ 

Vi = V,,^ -P^M, (7) 

CTn V Pn / 
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"(Pcxt+a2)2' 

Pn = \ {ol + Pext + B^j , (9) 

^=f^, (10) 
2pn 

gp = -Vp^, (11) 

V'-^-M-.^M/fc.], (12) 

-Bp = -Vp*, (13) 

* = [HBz,eq - S,ef)/fc,] . (14) 

In the above equations, Bp{x, y) = Bx{x, y)x + By{x, y)y is the planar magnetic field at 
the top surface of the sheet, Vn{x, y) — Vx{x, y)x + Vy{x, y)y is the velocity of the neutrals 
in the plane, and Vi(x,y) = Vi^x{x,y)x + Vi^y{x,y)y is the corresponding velocity of the 
ions. The operator Vp = x d / dx + y d / dy is the gradient in the planar directions within 
the sheet. The quantities tl){x,y) and '^{x,y) are the scalar gravitational and magnetic 
potentials in the plane of the sheet, derived in the limit that the sheet is infinitesimally 
thin. The vertical wavenumber fc^ = (fc^ + fc^)^/^ is a function of wavenumbers and 
ky in the plane of the sheet, and the operators T and represent the forward and 
inverse Fourier transforms, respectively, which we calculate numerically using an FFT 
technique. Terms of order 0{VpZ) in i^M, the magnetic force per unit area, are not 
written down for the sake of brevity, b ut are included i n the numerical code; their exact 



form is given in Sections 2.2 and 2.3 of ICiolek fc Basul ()2006l ). All terms proportional to 
VpZ are generally very small. 

The Eq. ([7]) above makes use of relations for the neutral-ion collision time Tni and the 
ion density rii, as described in BCW: 

ni\ + TO„ 1 

Tni = 1.4 ^ , , , (15) 

mi ni(o-w)iH2 

rii^/Cn^'. (16) 

Furthermore, Eq. ([7]) uses the normalized initial mass density (in units of cr„^o/io) Pn.o ~ 
1(1-1- Poxt), where Pcxt is defined below. 

Our basic equations contain three dimensionless free parameters: /io = 27rG^/^crn.o /^rcf 
is the dimensionless (in units of the critical value for gravitational collapse) mass-to-flux 
ratio of the reference state; fni^o = ''ni.o/^o is the dimensionless neutral- ion collision time 
of the reference state, and expresses the effect of ambipolar diffusion; i^xt = 2Pcxt /'I'Gcrg q 
is the ratio of the external pressure acting on the sheet to the vertical self-gravitational 
stress of the reference state. 

Each numerical model is carried out within a square computational domain spanning 
the region ~L/2 < x < L/2 and —L/2 < y < L/2. Periodic boundary conditions are 
imposed in the x and y directions. We choose a domain size L = IGttLq for all models 
presented in this paper, so that it is four times wider than the wavelength of maximum 
growth rate for an i sothermal nonmagnetic and unpressured sheet, AT.m = 47r_Lo (see 
Ciolek fc Basu 2006t BCW). The computational domain has N zones in each direction. 
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with N = 128 unless stated otherwise. Some results utilizing greater N are presented in 
Section [Ql 

Gradients in our simulation box are approximated using second-order accurate central 
differencing. Advection of mass and m agnetic flux is prescribed by using the mono- 
tonic upwind scheme of Ivan Leen (|l977^ . These forms of spatial discretization convert 
the system of partial differential equations to a system of ordinary differential equa- 
tions (ODE's), with one ODE for each physical variable at each cell. Time-integration of 
this large coupled system of OD E's is performed by using an Adams-Bashforth-Moulton 
predictor-corrector subroutine ( Shampin^ . Il994l ). Numerical solution of Fourier trans- 
forms and inverse transforms, necessary to calculate the gravitational and magnetic po- 
tentials ip and ^' at each time step (see Eqs. [12] and [H]), is done by using fast Fourier 
transform techniques. Further details of our numerical technique are given by BCW. 

The initial conditions of our model include a "turbulent" velocity field added to our 
background state of uniform column density anfl and vertical magnetic field strength 
i^rcf- T he velocity fi eld is generated in Fourier space using the usual practice adop ted 
by e.g. Stone et al. ( 1998h for three-dimensional models and Li fc Nakamura ( 20041 ) for 
thin-sheet models. For a physical grid consisting of N zones in each (x, y) direction, the 
wavenumbers and ky in Fourier space consist of the discrete values kj — 27rj/L, where 
j — ~N/2, ...,N/2. For each pair of k^ and ky, we assign a Fourier velocity amplitude 
chosen from a Gaussian distribution and scaled by the power spectrum oc fc", where 
fc — (fc^ -|- kyY^'^ . The phase is chosen from a uniform random distribution in the interval 
[0, 2tt]. The resulting Fourier velocity field is then transformed back into physical space. 
The distributions of each of and Vy are chosen independently in this manner, and 
each is rescaled so that the rms amplitude of the field is equal to Va- For n = —4, the 
resulting velocity field has most of its power in a large-scale flow component. We have 
experimented with various values of n and find that the results are generally similar 
as long as n < 0. Distinct results are found in the case of flat spectrum perturbations 
(n — 0), which we present for comparison. Finally, we note that velocity fields generated 
in the above manner are compressive. Hence, we have also studied one model with n — —4 
but the Fourier space amplitudes chosen so that Vx and Vy satisfy Vp • Vn = 0. 

Therefore, our turbulent initial conditions introduce the additional dimensionless free 
parameters Va/cg and n, while our simulation box introduces the parameters L/AT,m and 
grid size N. 



3. Results 



Table 1 contains a summary of models in our parameter study. For each model, we list 
the values of the free parameters fiQ, f^ifi, Pcxt, the form of the turbulent power spectrum, 
the turbulent velocity amplitude Va, the magnetosound speed T4is,o in the initial state of 
the model, and the calculated duration of the simulation, irun- The initial magnetosound 
speed is calculated from the other parameters of the mod el, and its relation to Va can 
act as a useful diagnostic. Following Ciolek fc Basu ( 2006h we use 
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Table 1 

Models and Parameters 



Model 


Mo 




Pcxt 


Spectrum 




Vms,o/cs 


trun/to 


1 


0.5 


0.0 


0.1 


k-^ 


2.0 


2.9 


> 5000 


2 


0.5 


0.2 


0.1 




4.0 


2.9 


0.8 


3 


0.5 


0.2 


0.1 


k-^ 


3.0 


2.9 


30 


4 


0.5 


0.2 


0.1 


k-^ 


2.0 


2.9 


31 


5 


0.5 


0.2 


0.1 


k-^ 


1.0 


2.9 


50 


6 


0.5 


0.2 


0.1 


k-^ 


0.5 


2.9 


58 


7 


0.5 


0.2 


10.0 


k-^ 


2.0 


1.0 


8.5 


8 


0.5 


0.2 


0.1 


fc-" (divO) 


2.0 


2.9 


23 


9 


0.5 


0.2 


0.1 




2.0 


2.9 


56 


10 


0.5 


0.1 


0.1 




2.0 


2.9 


92 


11 


0.5 


0.4 


0.1 




2.0 


2.9 


8.1 


12 


0.5 


1.0 


0.1 




2.0 


2.9 


1.8 


13 


1.0 


0.2 


0.1 




2.0 


1.7 


1.8 


14 


1.0 


0.2 


0.1 




1.0 


1.7 


4.3 


15 


1.0 


0.2 


0.1 


k-^ 


0.5 


1.7 


12 


16 


1.0 


0.2 


0.1 




2.0 


1.7 


33 


17 


2.0 


0.2 


0.1 




2.0 


1.2 


1.3 



where o ^/(l + Poxt) is the square of the normahzed (to Cs) initial Alfven speed, 

and = (1 + 3-Pcxt)/(l + ^cxt)^ is the square of a normahzed initial effective sound 
speed that includes the effect of external pressure, and follows from Eq. ((8]). Model 7 
has Cgg because for Poxt — 10, the large external pressure contributes significant 

opposition to the rest oring force of internal pressure in the initial state (see discussion in 
Ciolek fc Bas"ul . l2006l) . 



Our simulations are terminated as soon as cTn.max > 10 <Jnfl^ corresponding to runaway 
gravitational collapse of the first core. For models with Pcxt = 0.1, this also corresponds 
to a volume density enhancement by a factor « 100. We have verified with high resolution 
runs to greater values of crn.max/cn.o that the collapse does indeed continue. Collapsing 
regions are also invariably gravity-dominated, having locally supercritical mass-to-fiux 
ratios and net accelerations that point toward the density peak. This includes cases of 
prompt collapse, i.e. when collapse occurs in localized regions due to strong shocks as- 
sociated with the turbulent flow flcld, in a time t-c^-n less than 2 to- The values of irun do 
vary somewhat from one realization of the initial state to another, and in many cases 
represent an average value from many simulations. We have run each model at least 
five separate times, and some have been run significantly many more times, as described 
below. Model 1, which evolves under flux-frozen conditions (fni^o = 0), cannot be termi- 
nated in the usual manner. Since the initial mass-to-flux ratio is also subcritical for this 
model {^Q — 0.5), gravitational runaway collapse is not possible unless there is signifi- 
cant numerical magnetic diffusion. That simulation ran past t = 5000^0 without runaway 
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collapse or any notable artificial flux dissipation, thus providing an excellent verification 
of the accuracy of our code. While some models undergo many oscillations before even- 
tual runaway collapse of density peaks, others undergo prompt collapse. Although the 
models that undergo prompt collapse may be considered to be artificially forced into 
collapse by large-scale flows in the initial conditions, we nevertheless present them here 
as interesting limiting cases. Models 4, 13, and 17 constitute a special set of models with 
our standard neutral-ion couphng parameter fni,o = 0.2, external pressure Pext = 0.1, 
turbulent velocity amplitude Va = 2.0 Cg, but varying initial mass-to- flux ratio parameter 
fj,o = 0.5, 1.0, 2.0, respectively. Model 4 is run 15 times, and models 13 and 17 are run 25 
times (with unique random realizations of the initial velocity field), in order to compile 
statistics on the core mass distributions using the techniques described by BCW. Model 
4 can in some sense be considered our "standard" model since we are most interested 
in the acceleration of collapse in subcritical clouds due to nonlinear supersonic velocity 
perturbations. Model 8 is similar to model 4 but has the divergence-free initial velocity 
field. Model 3 is on the brink of either prompt collapse or a longer-term evolution lead- 
ing to runaway collapse, and can sometimes go into prompt collapse (see Section 4). For 
this reason, we ran the model over 15 times in order to yield a irun — 30 tg that is a 
characteristic value for all but a few runs that do go into prompt collapse. 

3.1. Global Properties 

Fig. [T] shows images of column density overlaid with column density contours and 
velocity vectors, for realizations of models 4 (top left), 13 (top right), and 17 (bottom 
left). Each snapshot is at the end of the simulation, when crn.max/cn.o = 10, but occurs 
at a different time irun, as indicated in Table 1. The maximum speeds in the simulation 
region are quite different at the end of the three simulations even though all three start 
with perturbations characterized by Va = 2.0 Cg. Therefore, the velocity vector plots each 
have a different normalization, with the horizontal or vertical distance between footpoints 
of vectors corresponding to 1.0 Cg, 2.0 Cg, and 3.0 Cg for the three models with /io — 0.5, 1.0, 
and 2.0 respectively. To understand the difference in maximum speeds, it is important 
to understand the different course of evolution in each model. 

The model 4, with /iq = 0.5, has a strong enough magnetic field that the initial 
compression driven by the large-scale flow of the nonlinear velocity field does not lead 
to prompt collapse in any region. The magnetic field causes a rebound after the initial 
compression. The densest regions never reexpand fully to the initial background density, 
and instead undergo oscillations in density until continuing ambipolar diffusion leads to 
the creation of regions of supercritical mass-to-flux ratio. These regions then undergo 
runaway collapse. For model 4, this occurs at a representative t — 31to, meaning that 
there is sufficient time for the initial velocity field to have damped significantly, since we 
do not replenish turbulent energy in these simulations. This is why the velocity amplitude 
is much smaller at the end of the simulation than in the other two runs. However, the 
maximum value is still supersonic (3.2 Cg), and there are strong systematic flow fields 
in the simulation. In contrast, when starting with small-amplitude initial perturbations 
(BCW), the maximum infall speeds are subsonic. In that case, runaway collapse also 
occurred much later, at t — 204 io- The color table and column density contours for 
model 4, when compared to those for models 13 and 17, reveal that the gas is not as 
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Fig. 1. Image and contours of column density <Jn{x, y)/(Tn,o, and velocity vectors of neutrals, for three dif- 
ferent models at the time that crn,max/o"n,o = 10- All models have g = 0.2, Pext = 0.1, and Va = 2.0 Cg. 
Top left: model 4 (/xq = 0.5). Top right: model 13 {^0 = 1.0). Bottom left: model 17 (/xo = 2.0). The color 
table is applied to the logarithm of the column density and the contour lines represent values of a"n/crn,o 
spaced in multiplicative increments of 2^/^, having the values [0.7,1.0,1.4,2,2.8,4.0,...]. The horizontal or 
vertical distance between the footpoints of velocity vectors corresponds to a speed 1.0 Cs for the /iq = 0.5 
model, 2.0 Cs for the /iq = 1.0 model, and 3.0 Cs for the /iq = 2.0 model. We use the normalized spatial 
coordinates x' = x/A-p ^ and y' = y/XT,im where Ax,m = 47r Lq is the wavelength of maximum growth 
rate from linear perturbation theory, in the nonmagnetic limit with Pext = 0. 



compressed and filamentary as in those cases, due to the rebound from the initial extreme 
compressions. 

The images, contours, and velocity vectors for models 13 (/io = 1.0) and 17 {fiQ = 2.0) 
reveal that the initial strong compression leads to immediate runaway collapse within 
highly compressed filaments. The velocity field has highly ordered supersonic compres- 
sive infall motions. The runaway collapse is occurring at times t = 1.8 io Sind t = 1.3 
respectively, essentially as soon as the initial flow creates a large-scale compression. Al- 
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Fig. 2. Image and contours of fi{x,y), the mass-to-flux ratio in units of the critical value for collapse. 
Regions with fi > 1 are displayed with a color table, while regions with /i < I are black. The contour 
lines are spaced in additive increments of 0.1. Left: Final snapshot of model 4 [fig = 0.5). Right: Final 
snapshot of model 13 (^o = 1.0). 

though kinetic energy is efficiently dissipated behind the shock fronts, there is hardly 
enough time for a large reduction of the global kinetic energy. Therefore, the maximum 
speeds at the end of the simulation (6.0 Cg for model 13 and 7.2 Cg for model 17) are quite 
similar to the initial maximum values. 

Fig. [2] shows images of the mass-to-flux ratio at the end of the simulations of models 
4 and 13. The end states have a combination of subcritical and supercritical regions. 
Subcritical regions are shown in black, and a color table is applied to the supercritical 
regions on both panels. The left panel illustrates that the supercritical regions of the 
initially significantly subcritical (fiQ — 0.5) cloud are created within the filamentary 
regions generated by the large-scale compressions. In contrast, the initially critical (/io = 
1.0) model has widespread supercritical regions generated by the small-scale modes of 
turbulence, as well as the most supercritical regions in the compressed layers. The former 
effect of widespread patches of mildly supercritical gas is possible due to the marginal 
nature of the critical (/io = 1-0) initial state. Physically, we might expect the cloud with 
[iQ = 1.0 to lead to a cluster of stars soon after the runaway collapse of the first cores. 
On the other hand, the significantly subcritical cloud with /io — 0.5 would show only 
isolated star formation in the compressed layers and have to wait a much longer time 
before ambipolar diffusion leads to clustered star formation in the remainder of the cloud. 

Fig.[3]shows the end states of models 4, 13, and 17 in different realizations (i.e. starting 
with a different but statistically equivalent initial velocity field) than shown in Fig. [TJ A 
color table shows the column density of the final state, and magnetic field lines above the 
sheet are also illustrated. These lines are generated in the manner described in BCW. 
The image of sheet surface and field lines above are viewed from an angle of 10° from 
the sheet normal direction. Animations of the evolution of the sheet surface density, with 
field lines appearing in the last frame, are available online. Clearly, the models which 
suffer prompt collapse (/ip — 1-0 and /ip — 2.0) show the most curvature of field lines 
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Fig. 3. Imago of gas column density crn{x, y)/a-n,o and superposed magnetic field lines for realizations of 
models 4, 13, and 17, with fig = 0.5 (top left), fio = 1.0 (top right), and = 2.0 (bottom left). All models 
have initial velocity amplitude Va = 2.0 Cg. These are two-dimensional projections of three-dimensional 
images containing a sheet with a column density image and magnetic field lines extending above the 
sheet to a distance about half the box width. The image is seen from a viewing angle of about 10° 
relative to the sheet normal direction. Animations of the evolution of the column density are available 
online. The field lines appear in the last frame of the animation. 

since the field is dragged inward by the strong initial compression wave. In contrast, 
the cloud with /ig — 0.5 (top left) undergoes a rebound and several oscillations before 
runaway collapse can occur. This allows the magnetic field to straighten out again. The 
ultimate collapse of the first core is due to ambipolar drift of neutrals past field lines, 
so the field is not significantly distorted by this process. However, a legacy of the initial 
compression is that the mass-to-flux ratio is no longer spatially uniform, and significant 
column density structure exists even if the magnetic field is not very distorted. 

The relative amounts of field line curvature in the cloud and within dense cores are 
quantified by 6* = tan"i(|Bp|/B^,oq), where \Bp\ = {Bl+B^y/"^ is the magnitude of the 
planar magnetic field at any location on the sheet-like cloud. Hence, 6 is the angle that a 
field line makes with the vertical direction at any location at the top or bottom surface of 
the sheet. To quantify the differences in field line bending from subcritical to transcritical 
to supercritical clouds, we note that representative realizations of models 4, 13, and 17, 
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x' x' 

Fig. 4. Top; Column density and velocity vectors as in Fig. [T] but for models 9 and 16, which have initial 
nonlinear velocity field with Va = 2.0 Cs and flat power spectrum (k^). The horizontal or vertical distance 
between the origins of velocity vectors corresponds to a speed 0.5 Cg. Bottom: Images of mass-to-flux ratio, 
as in Fig. [2] but for models 9 and 16. 



with /io = (0.5, 1.0, 2.0), have average values (?av = (4.2°, 18°, 23°), and maximum values 
(probing the most evolved core in each simulation) 0max = (25°, 65°, 65°). Of these, only 
the model with /ip — 0.5 shows similar values of as a corresponding model with initial 
small-amplitude perturbations. For models with small-amplitude initial perturbations, 
BCW found that = (0.5,1.0,2.0) yields representative values e^av = (1.7°, 8.3°, 18°), 
and 0„,ax = (20°,30°,46°). 

The evolution of models with flat spectrum (u^ oc kP) initial perturbations is distinct 
from the cases with negative exponent, so we present the results from models 9 and 16 in 
Fig. m In these models, Va = 2.0 Cg as in model 4, but the large-scale flow does not domi- 
nate the initial condition. Therefore, the small-scale modes contribute more significantly 
to enhance ambipolar diffusion. This enhancement of ambipolar diffusion due to small- 
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scale i rreg ularities i s simil ar to the mechanism studied analytically by Fatuzzo fc Adam"^ 
(|2002t ) and lZweibell (|2002f) . We study only models with /ig = 0.5 and /xq = 1.0 in order to 
focus on the enhanced ambipolar diffusion. The values of irun for these models are 56 to 
and 33 to respectively. These are significantly longer time scales than in the corresponding 
models with oc k~^. However, they are much shorter than in models with the same 
background state and linear initial perturbations (studied by BCW), in which 
is 204 to and 121 to for /io ~ 0.5 and — 1.0, respectively. 

Fig. |4] shows the column density and mass-to-flux ratio at the end of simulations of 
model 9 and 16. The input turbulence acts to increase the rate of ambipolar diffusion, 
but the turbulence also decays away. By the time of runaway collapse, the cloud structure 
and kinematics more closely resembles the case of small-amplitude initial perturbations 
than the case of nonlinear-flow-induced fragmentation (models 4 and 13). Representative 
values of the maximum speed Wmax of neutrals at the time t — tiun of runaway collapse 
are 0.7 Cg for model 9 and 0.8 Cg for model 16. Both values are closer to the values 0.4 Cs 
and 0.7 Cg when starting with small-amplitude initial perturbations (BCW) than for the 
cases of nonlinear-flow-induced fragmentation, in which case Wmax — 3.2 Cg and 6.0 Cg, 
respectively. The bottom panels show the mass-to-flux ratio at the end of simulations of 
model 9 and 16. The subcritical model 9 has only isolated pockets of supercritical cores, 
as well as emerging cores which still have subcritical but enhanced mass-to-flux ratio. 
The image is similar to the corresponding image when starting with small-amplitude 
perturbations (Fig. 9 of BCW), but the cores are not circular in shape. The fragmentation 
scale also seems related to that of the small-amplitude perturbation model, although 
many fragments are just beginning to emerge and may take a much longer time to develop 
fully. The corresponding image for model 16 shows that the initially critical (/zq = 1.0) 
state leads to many regions of supercritical mass-to-flux ratio. The emerging fragment 
pattern has less resemblance to the corresponding case that starts from small amplitude 
perturbations (Fig. 9 of BCW) than for /io = 0.5. Nevertheless, a fragmentation scale is 
more apparent than for the case of nonlinear-flow- induced fragmentation shown in Fig. [21 
The relatively low amount of remaining turbulence at the end of both simulations means 
that new cores will develop on the non-accelerated ambipolar diffusion time scale, leading 
to a significant age spread of star formation. In the nonlinear initial condition models, 
with k^* or /c" spectrum, we may identify the initial cores with an early phase of star 
formation that is induced by turbulence, and the emerging cores with a later phase of 
star formation that can grow from lingering small-amplitude perturbations. In this way, 
our mode l could be in qualitative agreement with the empirical scenario advocated by 
IPalla fc Stahleij (|2000h . 



3.2. Time Evolution 

Fig. E] shows the time evolution of maximum column density for four models and 
of the maximum mass-to-flux ratio for three models, all having fj,o = 0.5. It clearly 
shows that for subcritical clouds: (1) fragmentation by runaway collapse does not occur 
under conditions of flux-freezing (dash-dotted line; the simulation actually runs past 
t = 5000 to without runaway collapse); (2) small-amplitude (linear) perturbations result 
in the classical quasistatic evolution requiring a time tiun ~ 200 to; (3) flat spectrum (fc°) 
nonlinear perturbations that result in a collapse time that is shorter by a factor w 4; (4) 
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Fig. 5. Time evolution of maximum values of surface density and mass-to-flux ratio for various models 
with fiQ = 0.5. The solid lines show the evolution of the maximum value of surface density in the simu- 
lation, (Tn,max/cn,Oi versus time t/to. The dashed lines show the evolution of the maximum mass-to-flux 
ratio in the simulation, ^max. This is shown for models 4, and 9, which have = 0.5 and same values 
of fjii and Pexti but different power spectra of turbulent initial perturbations, k~'^ and fc" respectively. 
Two other models are also shown for comparison. One has the same parameters but linear initial pertur- 
bations, corresponding to model 1 of BCW. Furthermore, the dash-dotted line shows the evolution (up 
to time t Ri 200 to only) of (Tn^max/crnjO for the flux-frozen (rni^o = 0) model 1, which never undergoes 
runaway collapse. 

power-law {k~^) spectrum nonlinear perturbations that result in a rapid collapse that 
is shorter than the linear case by a factor « 7. Of course, the exact values of t,-un will 
depend on Va and other parameters such as fni^o and Pcxt- For the fourth case above, 
it may depend on the box size L as well, since that sets the scale of th e largest mode 
in the simulation. This figure corresponds to Fig. 1 of Kudoh fc BasrJ (|2008,) . which 



shows results for some three-dimensional models. In their figure, the volume density 
Pn and plasma f3 (counterparts to an and fi in our thin-sheet model) undergo some 
variations due to vertical oscillations of the cloud, before eventually increasing rapidly. 
Our model follows the integrated quantities through the layer and therefore does not 
include the effect of vertical oscillations. Nevertheless, the timescale of evolution and 
eventual runaway collapse are in good agreement where comparisons can be made. Our 
thin-sheet model allows a broader parameter study than currently possible using three- 
dimensional simulations. 

Fig. [6] shows the time evolution of maximum column density and maximum mass-to- 
flux ratio for three models with fiQ ~ 1.0. The model with small amplitude (linear) initial 
perturbations corresponds to model 3 of BCW. The other ones are model 13 and model 
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Fig. 6. Time evolution of maximum values of surface density and mass-to-flux ratio for various models 
with ^0 = 1-0. The solid and dashed lines have the same meaning as in Fig. [5] Shown are results from 
models 13 and 16, which have differing power spectra of turbulent initial perturbations. Also shown is a 
model with the same parameters but linear initial perturbations, corresponding to model 3 of BCW. 

16 of this paper. The same three quahtatively distinct evolutionary modes occur as for 
the case of /Zq = 0.5. A notable difference is that collapse occurs immediately during the 
first compression in the case of nonlinear- flow- induced fragmentation, at t^^m = 1.8 io- 
There is no rebound from the first compression as occurs when /zq = 0.5. 

Fig. [7] illustrates the effect of varying levels of ionization on the fate of nonlinear- 
flow-induced fragmentation. This is represented by differing values of fni^O: with fni^o = 
corresponding to flux-freezing and differing values of fni,o corresponding to different initial 
ionization fractions oc Tj^q (see Appendix). The standard model with fni,o = 0.2 
corre sponds to a canonical ionization fraction Xi^o ~ 10~^(nn.o/10^ cm"'^)"^/^ ( Tielend . 



20051 ). Our results show that ^run may indeed vary significantly due to variations in 
Xifi, easily spanning the range of 10^ yr to 10^ yr for typical values of units and input 
parameters. Clearly, a definitive understanding of the influence of magnetic fields awaits 
further insight into ionization levels in molecular clouds. 



3.3. Core Properties 

Two important observed properties of dense cores are the kinematics of infall mo- 
tions, and the distribution of core masses. The latter suffers from some ambiguity due 
to the different possible definitions of a "core" , or more specifically, how to define a core 
"boundary" . Here we describe the most basic features of our simulated cores. 
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Fig. 7. Effect of differing levels of magnetic coupling on time evolution of the maximum value of surface 
density. The solid lines show results for models 4, 10, and 11. All models have no = 0.5, Pext = 0.1, and 
turbulent initial perturbations with power spectrum oc fc~*. The dashed line shows the evolution of the 
flux-frozen model 1 o = 0), which does not undergo runaway collapse. 

Fig. [5] shows the velocity profiles (using dashed lines) in the vicinity of cores that 
are obtained in simulations of models 4, 13, and 17, with fiQ = 0.5,1.0, and 2.0, re- 
spectively. They are measured along a line in the a;-direction that passes through the 
column density peak. Also shown (in dotted lines) for each value of fiQ is a profile of 
x-velocity through a core that is formed in a simulation with small-amplitude (linear) 
perturbations (BCW), but otherwise the same parameters as the other model shown in 
the same panel. The horizontal solid lines in each panel pass through the "midplane" 
of the velocity profiles, and allow one to read off the systematic x- velocity of the core. 
For the cases of initial linear perturbations, the increasing sequence of fiQ leads to ever 
increasing maximum infall speeds, from about half the sound speed up to mildly super- 
sonic values. There is also evidence of infall speeds increasing towards the core centers, 
due to gravitational acceleration. For the cases of initial nonlinear perturbations (these 
models all have a spectrum), the sequence of increasing leads to greater relative 
infall speeds onto the cores. These motions are supersonic in all cases, and constitute an 
important observationally-testable consequence of nonlinear-flow-induced fragmentation. 
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Fig. 8. Velocity component of neutrals Vx, normalized to Cs, along a line parallel to the x-axis that passes 
through the center of a core, for models with various values of fig and differing initial perturbations. The 
horizontal coordinate x' = {x — Xc)/Lo, where Xc is the location of the core center in each case. The 
dashed lines show the profile of Vx for cores generated in models 4, 13, and 17, from left to right panels. 
They are characterized by initial nonlinear perturbations (NLP). The dotted lines show for comparison 
the profiles through cores in models 1, 3, and 5 of BCW, i.e. models with the same parameters but initial 
linear perturbations (LP). The horizontal solid lines denote the systematic speeds of the cores in the 
X— direction. Note the largest systematic speed in the model with fiQ = 0.5 and initial NLP. 

The models with greater fiQ have greater infaU speeds because they undergo collapse dur- 
ing the first compression, with most of the initial input turbulent energy still intact, i.e. 
there has not been much time for turbulent decay. Also note that there is no evidence 
for an accelerating flow in these cases, which would be a signature of gravitationally- 
driven motions. Thus, these models demonstrate flow-driven core formation, rather than 
gravitationally-driven core formation. For the model with /io = 2.0 there is essentially 
no systematic core speed, since the collapse occurs very quickly at the intersection of 
two colliding flows. At the other limit of a significantly subcritical cloud (/io = 0.5), the 
initial compression is followed by a rebound due to the strong magnetic restoring forces. 
The core forms later within the region of high density that is undergoing oscillatory 
motions. The systematic speed of the core relative to the simulation box is supersonic 
(about twice the sound speed in this model), although the relative speed of infall onto 
the core is subsonic or transonic. The large systematic core speeds for subcritical clouds 
constitute another observationally-testable consequence of nonlinear-flow-induced frag- 
mentation. The case of fia — 1.0 is intermediate in features between the fxo = 0.5 and 
fio = 2.0 models, but is actually closer to the fio = 2.0 model since the collapse occurs 
very quickly, with — 1.8 to- This is very close to the value trun = 1.3 to for the 
/io = 2.0 model. 

Fig. [S] shows the histograms of core masses, defined as masses enclosed within regions 
that have tTn/cnjO ^ 2 surrounding a column density peak. These are measured at the 
end of each simulation, for 15 separate realizations of model 4, and 25 each of models 
13 and 17. Simulations of model 4 (/io = 0.5) and model 13 (/io — 1.0) produce an 
average of five identifiable cores per simulation, while model 17 (/^o = 2.0) produces an 
average of ten cores per simulation. For details about our thresholding technique used to 
obtain core masses, see BCW. The histograms reveal that for any fixed value of /io, the 
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Fig. 9. Histograms of masses contained within regions with fTn/un^o ^ 2, measured at the end of 
simulations with parameters of models 4, 13, and 17. Specifically, they are distinguished by values of 
/iQ = 0.5, 1.0, 2.0 as labeled. Each figure is the result of a compilation of results of many simulations. 
The bin width is 0.1. 



distribution of core masses is much broader than the corresponding histogram of masses 
for fixed /zq and initial small-amphtude (hnear) perturbations. See Fig. 8 of BCW for 
the latter, which show a very sharp descent beyond the preferred mass scale. There are 
many more high-mass cores that are formed in these models. However, an examination of 
Figs. [Hand [3] reveals that many of the cores have very elongated and irregular shapes, so 
that they may yet break up into multiple fragments. In BCW, we proposed that a broad 
CMF may be caused by a distribution of initial mass-to-flux ratio values within a cloud. 
That remains an alternative scenario to that of "turbulent fragmentation" explored here 



and in several previous publi cations (jPadoan et al.l . 119971 : iKlessenl . 1200 It iGammie et al 
20031 : iTillev fc PudritzL l2007h . 



3.4. Turbulent Dissipation 



The rate of dissipation of turbule nt energy has been studied extensively in a series of 



three-dimensional simulat ions (e.g. Stone et all . 1 19981: Mac Low et al 



1999; 



Ostriker et al.l.l200lh. See also th e reviews bv lMac Low fc Klessen 



1998 



mm 



Mac Low 



Elmegreen fc Scalol 



(|2004[ ). and lMcKee fc Ostriken (j2007[ ). In this Section, our goal is to briefly present some 
information about the turbulent decay in our simulations. These may be of interest be- 
cause our simulations differ in their use of the thin-sheet approximation and the use of 
high-order adaptive time-stepping that is part of our implementation of the method of 
lines technique (see BCW). Nevertheless, we do also obtain relatively rapid turbulent 
dissipation in most models, as presented in the various figures in this Section. 

We present results for the decay of kinetic energy in our simulation box. However, 
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Fig. 10. Kinetic energy versus time for models with various values of no and/or Va (normalized to Cs). 
Each of these models is run with A'^ = 256. 

the total energy in the simulation box is not conserved, due to radiative losses implied 
by our isothermal assumption, and also due to work done by the external pressure and 
magnetic forces associated with the field components and By at the cloud's top and 
bottom surfaces. Nevertheless, the amount of turbulent kinetic energy present in a cloud 
has important observational implications, so we illustrate its evolution here in several 
figures. We also use the kinetic energy evolution as a means of exploring the effect of 
numerical resolution in our simulations. 

Fig. [TUl shows the evolution of kinetic energy i?kin (defined as the sum of |o'n(u^ + Vy) 
over all cells of a simulation), normalized to initial values £'kin,0: for models 4, 6, and 13. 
Note that the values of i?kin,o differ from model to model. The model with fiQ — 1.0 does 
not have much chance to lose kinetic energy because collapse occurs right away, during 
the first turbulent compression. For models with fio = 0.5, there is a rebound from the 
initial compression, and this is indicated by the oscillations of £'kin- Furthermore, there 
is an overall systematic decay of E'kin so that it is significantly reduced in one sound 
crossing time of the initial half-thickness of the cloud, tc = Zq/cs — 2Lq/cs = 2tQ, where 
we have used Eq. (30) of BCW to relate Zq to Lq. The decaying oscillations of i?kin 
are consistent with the qualitative picture obtained from the animation of model 4 that 
accompanies Fig. [31 Some of our realizations do show an increase in i?kin during the last 
stage of evolution, due to the conversion of gravitational energy into kinetic energy of 
systematic infall onto one or more cores. 

Fig. [m shows the effect of resolution on the decay of kinetic energy. Our standard 
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simulations have N = 128, and numerical experiments with = 256 and = 512 
demonstrate that while some more kinetic energy is retained in those cases, the overall 
pattern of decay and oscillations of -Ekin is maintained. We note that each simulation has 
a unique random but statistically equivalent initial state. 

Fig. [T2l reveals the additional effect on the kinetic energy evolution of two interesting 
limits. In one case, we perform the numerical experiment of starting with a divergence- 
free (non-compressive) initial velocity field, e ven though turbulence in the interstellar 
medium is thought to be highly compressive (jMcKee fc Ostrikeii 120071 ). For this case, 
a plot of £^kin for the evolution of model 8 (dashed line) reveals that the kinetic en- 
ergy still decays, but that the cloud does not undergo large-scale oscillations during the 
process. These oscillations occur in the compressible case due to the restoring force of 
the magnetic field when compressed into filamentary structures. This does not occur in 
the incompressible model in a globally coherent manner, although locally compressive 
motions are generated during the evolution and turbulent decay does occur rapidly, as 
in the other models. The dash-dotted line reveals the interesting evolution in the case 
of flux- freezing (model 1; fni^o = 0)- Here the initially compressive velocity field leads 
to an initial rapid decay of turbulence through shocks, but large-scale oscillatory modes 
remain in the simulation box for indefinite periods of time. These modes have a root 
mean squared velocity amplitude « 2 Cg and contain roughly half the initial input energy. 
Why do these modes not decay away? The lack of ambipolar diffusion means there is 
no dissipation of modes in which the restoring force is due to the magnetic field. The 
k^^ spectrum means that the largest modes dominate, and these also suffer negligible 
numerical dissipation in our scheme. The restoring force that drives the waves is provided 
largely by the magnetic tension associated with the magnetic field external to the sheet. 
While the case of a thin sheet may not be generalizable to three dimensional molecu- 
lar clouds, we feel this result is an important pointer to processes that may in fact be 
occurring in real clouds. That is, the external magnetic field, anchored in the Galactic 
interstellar medium , may allow the outer parts of clouds (effectively fiux-frozen due to 



• parts c 

UV ionization - see ICiolek &: Mouschovia to maintain long-lived oscillations that 



ar e then identified observati onally as "turbulence". This idea has long been advocated 
bv iMouschoviad ( 1975L 1983) ■ We note that this result could not be obtained in periodic 
box simulations that contain no effect of an external medium, and leave a more thorough 
assessment of this effect to a forthcoming paper. 



4. Discussion 



We have performed a parameter study of fragmentation of a dense sheet aided by the 
presence of initial nonlinear velocity perturbations. In most models, the power spectrum 
of fluctuations is oc fc^*, so that the initial conditions impose primarily a large-scale 
flow to the system. We have also studied the case of nonlinear perturbations with power 
spectrum cx in which the small-scale fluctuations play a bigger role. Of the two modes, 
the latter is more similar to gravitational fragmentation arising from small-amplitude 
perturbations, as studied extensively in our previous paper (BCW). The main difference 
is an accelerated time scale for core formation. This is particularly apparent for the 
cases with subcritical initial ma ss-to-flux ratio, in which case the nonlin ear fluctuations 
enhance ambipolar diffusion (see Fatuzzo &: Adanii . 20021 [Zweibel [200^ . For the case of 
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Fig. 11. Kinetic energy versus time for model 4 parameters but varying resolution. 



nonlinear-flow-induced fragmentation, originally studied bv iLi fc Nakamural ( 2004( l and 
Nakamura fc Lil (j2005r ). we find that the induced structures are highly filamentary and go 
into direct collapse for supercritical clouds. For critical, and more so for subcritical clouds, 
the initial compression may be followed by a rebound and oscillations which eventually 
lead to runaway collapse in dense pockets where enhanced ambipolar diffusion has created 
supercritical conditions. What determines the outcome? For any given field strength, 
there is a threshold initial velocity amplitude Va above which prompt collapse will take 
place. An examination of model outcomes in Table 1 reveals that, for a fixed standard 
initial ionization fraction defined by fni^o = 0-2, prompt collapse takes place when Va > 
^MS,o (see models 2, 13, and 17). Indeed, the model 3, which has Va ~ Vms,0: is actually 
prone to go into prompt collapse (trun ~ ^o) in some realizations, but undergoes several 
oscillations before runaway collapse in most cases (with representative value irun = 30 ip)- 
We can say that significantly super- Alfvenic perturbations are associated with prompt 
collapse, for both subcritical and supercritical model clouds. This criterion does not 
apply to models with initial power spectrum oc since the kinetic energy does not 
get channeled toward a large-scale compression wave. It also does not apply to model 7 
(^ext = 10), since its low value of Vms,o is very specific to the external-pressure-dominated 
initial state, but not representative of the signal speed in the high-density regions that 
are subsequently generated. 

The highly filamentary structure of clouds in which prompt collapse takes place is a 
source of concern when comparing with maps of observed molecular clouds. This was 
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Fig. 12. Kinetic energy versus time for models 1 (dash-dotted line), 4 (solid line), and 8 (dashed line). 
Each of these models has fig = 0.5 and Va = 2cb and is run with A'^ = 128. They differ in their values of 
fni and in whether or not the initial velocity field is divergence-free. 



noted bv lLi &: Nakamura ( 2004 ). who commented that clouds with weak magnetic field 
and supersonic turbulence as modeled (i.e. fc~^ power spectrum, having most power on 
the largest scales) would appear too filamentary in comparison with observations. Our 
study extends this concern also to models with strong magnetic field if the turbulence 
is highly super- Alfvenic, since prompt collapse occurs in highly compressed filaments 
without a chance for them to rebound. Since the weak magnetic field cases also have 
by design a velocity amplitude that is super- Alfvenic, we can say that super- Alfvenic 
turbulence in all cases may have problems with excessive filamentarity. There is another 
problem with large amounts of turbulent forcing; the relative infall motions onto the 
cores are high l y supersonic (see Fig . Et, and at odds with observed core i nfall motions 



(Tafalla et al.. 1998: Williams et al 



s. 01), ana at odds witn observed core i mall motions 
D, fl999t iLee et all . 120011 ICaselh et aP . l2002f ). which 



are subsonic or at best transonic. Of course, both of these problems are set up artificially 
in our simulations through nonlinear forcing associat ed with the initial conditio ns. In 
other simulations of driven super- Alfvenic turbulence ( Padoan Sz Nordlund . 20021 ). such 
forcing continues at all times and the above features are always present. 

If the highly turbulent and/or super- Alfvenic models pose difficulties for dense core 
formation, then how does one account for the highly supersonic motions observed in 



molecular clouds fe.g. ISolomon et al.l . 119871 )? The answer is likely that they exist in the 



lowest density envelopes of the molecular clouds and therefore should not be input into 
local models of dense subregions, as we do in some cases here. Our super- Alfvenic models 
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in a periodic simulation box demonstrate a limiting case, and help establish that such 
models cannot be applied directly to explain observed star-forming regions. In a global 
scenario, the dense regions that form cores will be less turbulent than the larger low- 
density envelopes. The low-density regions can support highly turbulent motions while 
denser regions have lower velocity disp ersion, as demonstrated in 1 .5 dimensional global 
models of molecular cloud turbulence ( Kudoh fc 

Our Fig. [5] shows that the core mass distribution is relatively broad for any given 
value of /io for nonlinear-flow-induced fragmentation (fc~^ spectrum of velocity fluctua- 
tions). This repeats t he quali t ative findings of many earUer studies in thre e -dime nsions 
( Padoan et al. . 1997 : KlessenL 2001 : Gammie et al. . l200,l iTiUev fc Pudritj - lioorh . This 
scenario of turbulent fragmentation is a plausible mechanism to generate broad CMFs 
of the type observed. It remains an open question whether this kind of CMF is related 
to the IMF since the cores are often highly irregular in shape, and it is not clear that 
they will collapse monolithically. Alternative methods to generate broad IMFs or CM Fs 
are the global effect of competitive accreti on (Bonnell et al. . 20031 Bate et al. . 2003), a 
temporal spread of core accretion lifetimes ( Mversl . l2000l : Basu fc JoneZ ^ 2004( l. or a distri- 
bution of initial mass-to-flux ratios in a cloud (BCW). Future work by the astrophysical 
community may clarify the relative roles of these processes. 



5. Summary 

We have studied the effect of initial nonlinear velocity perturbations on the formation 
of dense cores in isothermal sheet-like layers that may be embedded within larger molec- 
ular cloud envelopes. Our simulation box is periodic in the lateral (x, y) directions and 
typically spans four nonmagnetic (Jeans) fragmentation scales in each of these directions. 
The initial input turbulent energy is allowed to decay freely. The simulations reveal a 
wide range of outcomes. We emphasize the following main results of the paper: 

(i) Time Evolution to Runaway. Subcritical model clouds can undergo accelerated am- 
bipolar diffusion in two different ways. For nonlinear initial velocity perturbations 
in which small-scale modes contain a large portion of the energy, the onset of run- 
away collapse occurs sooner by a factor « 4 in our typical models. For nonlinear 
perturbations with most energy on the largest scales (hereafter, nonlinear ffows), 
the runaway collapse can be sped up by a greater factor, w 7 for our typical models. 
Supercritical clouds undergo prompt collapse whenever nonlinear ffows are present. 
Subcritical model clouds may also be pushed into prompt collapse by nonlinear 
flows that are significantly super- Alfvenic. 

(ii) Morphology of Clouds. Supercritical model clouds whose evolution is initiated by 
nonlinear flows have a highly fllamentary structure. Subcritical clouds with initial 
nonlinear but trans- Alfvenic or sub- Alfvenic flows have a markedly less filamentary 
structure. In these cases, magnetic fields cause a rebound from the initial compres- 
sion, and several oscillations occur before the runaway collapse of the first cores. 
Subcritical clouds with initially super-Alfvenic nonlinear flows promptly develop 
highly filamentary structure with embedded collapsing cores. 

(iii) Velocity Profiles. Supercritical and transcritical model clouds which are driven into 
prompt collapse have highly supersonic infall speeds at the core boundaries, while 
subcritical model clouds typically have transonic or subsonic infall speeds (relative 
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to the velocity centroid) onto cores. In the subcritical cases, the cores can have larger 
systematic motions than in supercritical models, because the cores form within 
regions undergoing oscillatory motions. We believe that the large infall motions 
in the models with super-Alfvenic nonlinear flows may disqualify them as viable 
models for core formation, given current observational resiilts. 

(iv) Core Mass Distributions. Core formation initiated by nonlinear flows leads to 
broader core mass functions than found in earlier studies of fragmentation initi- 
ated by small-amplitude perturbations. This applies to models of any fixed initial 
mass-to-flux ratio /iq. However, the ultimate relation of such a core mass function 
to the stellar initial mass function is not settled due to the irregular shape of the 
fragments created by nonlinear flows. These fragments may in turn break up into 
multiple objects at a later stage. 

(v) Turbulent Decay. Supersonic initial velocity perturbations lead to an initially rapid 
decay of kinetic energy in all models, on a time scale similar to the sound crossing 
time across the half-thickness of the sheet. This rapid decay of turbulence is in 
agreement with a wide variety of previous results in the literature. However, sub- 
critical model clouds can undergo oscillations that reduce the decay rate of kinetic 
energy at later times. Furthermore, in the limit of excellent neutral-ion coupling 
(flux- freezing), as may be present in UV- ionized molecular cloud envelopes, large- 
scale wave modes may survive for very long times. 
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Appendix A. Typical Values of Units and Other Quantities 

Typical values of our units are 

1 /2 

Cs = 0.188 (^^^ kms-i, (A.l) 



10 k; V 

Lo = 7.02 X 10-^ [j^ji ^ , ) PC (A.3) 



*n,0 



= 1.45x10MJ^1 r^°" I AU, 



10 k; V ^n,o 

„„ = 9.19xlO->(j^)'(i2;i^^lMe, ,A.4) 



24 



Here, wc have used N^fl = o'n.o/'Ttn, where = 2.33 mpj is the mean molecular mass of 
a neutral particle for an H2 gas with a 10% He abundance by number. Furthermore, we 
may calculate the number density of the background state as 

n„,o = 2.31 X 10^ (HiS) ' (1 + Pext) cm-3. (A.6) 

The dimensional background reference magnetic field strength for a given model is simply 
Bref = Bq/hq. Finally, the ionization fraction (= ni/rin) in the cloud may be expressed 
as 
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